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Fast Inverse Dose Optimization (FIDO) for IMRT via 
Matrix Inversion without Negative Beamlet Intensities 
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Abstract 

A fast optimization' algorithm is very important for inverse 
planning of Intensity Modulated Radiation Therapy (IMRT), and 
for adaptive radiotherapy of the future. Conventional numerical 
search algorithms such as the conjugate gradient search, with 
positive beam weight constraints, generally require numerous 
iterations and may produce suboptimal dose results due to 
trapping in local minima. A direct solution of the inverse problem 
using conventional quadratic objective functions without positive 
beam constraints is more efficient but will result in unrealistic 
negative beam weights. We present here a direct solution of the 
inverse problem which does not yield unphysical negative beam 
weights. The objective function for the optimization of a large 
number of bearalets is reformulated such tiiat the optimization 
problem is reduced to a linear set of equations. The optimal set of 
intensities is found through a matrix inversion, and negative 
beamlet intensities are avoided without the need for externally 
imposed ad-hoc constraints. The method has been applied to a test 
phantom and to a few clinical radiotherapy problems. We achieve 
highly confonnal primary dose distributions with very short 
optimization times. Typical optimization times for a single 
anatomical slice (2D) using a Matlab matrix inversion routine in a 
single processor desktop computer, are: 0.03 sec. for 400 
beamlets; 0.25 sec. for 1,000 beamlets; 3,5 sec. for 2,000 beamlets 
and 8 sec for 3,000 beamlets. In conclusion, the new method 
provides a fast and robust technique to find a global minimum that 
yields excellent results for the inverse planning of IMRT. 

Keywords: Inverse planning, optimization, objective functioiL 
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1. Introduction 

Intensity Modulated Radiation Therapy (IMRT) is becoming a new standard technique for 
radiotherapy » Given the availability of online imaging tools (Mackie et al 1999, Groh et al 
2001, Jaffray et al 2002) and the improved conformal dose distributions obtained through 
IMRT and its dynamic delivery features, adaptive radiotherapy becomes an important factor to 
be considered. A fast and reliable optimization algorithm is crucial not only for designing good 
radiation treatment plans but also for the successful implementation of future interactive 
adaptive treatment techniques. Conventional optimization algorithms using numerical searches 
such as the conjugate gradient search (Smith and Lasdon 1992, Spirou and Chui 1998) with 
positive beam weight constraints usually require numerous iterations (i.e. long computation 
times) and may result in suboptimal plans due to trapping in local minima of the objective 
function. A direct solution of the inverse problem using conventional quadratic objective 
functions without imposing positive beam constraints would be computationally faster but will 
result in unrealistic negative beam weights. We present here a fast inverse dose optimization 
(FIDO) method for the direct solution of the inverse problem that avoids negative beamlet 
weighs and preserves efficiency. 

The most fundamental requirements of a radiation treatment optimization are: (i) dose is 
homogeneously deposited in the Planning Target Volume (PTV); («) the dose deposited in any 
Organ At Risk (OAR) does not exceed a threshold value and ideally should be zero; (iii) the 
dose deposited in the rest of the organs and tissue not included in the PTV and OARs, should 
be as low as is reasonably achievable (ALARA principle) to minimize the risk of secondary 
carcinogenesis- We will call these global whole-body regions "all the rest", abbreviated as 
ATR; (iv) the planned dose gradient crossing the PTV boundaries should be as high as 
possible, but it is subsequently degraded by practical considerations such as patient setups. 
Optimizations are pursued by the minimization of a positive-definite objective function. A 
successful optimization will yield a global minimum to this objective function in a short 
computation time with physically achievable beamlet intensities. 

This work reintroduces an approach to radiation treatment optimization that is very fast by 
solving a linear system of equations and yields a global minimum of the objective function 
without the use of a search routine. The possibility of such a direct optimization scheme has 
been known for decades but it has been "impossible** to implement rigorously using 
conventional quadratic objective functions because optimized results are only achievable with 
unphysical negative beamlet intensities (Webb 2003). If an ad-hoc condition requiring tiie 
beamlet intensities to be positive is introduced (e.g. force negative values to be zero), the linear 
solution method yields a dose distribution with artefacts and non-uniformities. Alternatively, 
one may pursue a minimization of the objective function by a direct search over all possible 
positive beam intensities, but this is a very time-consuming "trail and error" approach. 
Practically, such approaches become prohibitive if repeat optimizations would be desirable but 
avoided; the dose results are sometimes deemed "good enough" prematurely. 

In this work the negative-intensities problem is considered within the objective fimction 
and not as an externally imposed ad-hoc requirement. In this way the objective function 
remains quadratic and the single global minimum is obtained by efficient matrix inversion. 
Only primary KERMA has been included in the calculations reported here; however an 
extension to full 3D dose optimization is straightforward. 

LI Convjsntional objective Junction 

For simplicity of the formulation we will consider a single PTV, a single OAR and the ATR 
region. A typical objective function 6 satisfying the optimization conditions stated above is of 
the form: 



ATR 



(1) 



where the pfc are importance coefficients and the objectivity terms are: 
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0^=1.{ Z MW-rf'"' . - (2a) 

OoAR = Z Z ^i^M > (2b) 

and <^.»= Z Z ♦^/''rW . (2c) 

where is the weight of beamlet /, is the dose deposited at destination point x by beamlet i 
and df^ is the dose prescribed to the PTV. For the moment and for the sake of clarity, the 
importance coefficients in the objective function 6 (eq. 1) are constants. Later in tiie paper 
(Section 3.1) we will generalize this technique to the case of position-dependent importance 
coefJScients (Shou and Xing 2004). 

The objective function O (eq. 1) is positive and has an absolute minimum equal to zero 
when each of the terms is zero. This means a dose equal to tf^^ at each point inside the PTV 
and a zero dose elsewhere. The main reason for the traditional appearance of negative weights 
upon optimization of the objective function O is the fact that we inq>ose two conflicting 
conditions: we require 6^^^^ = 6^^ = 0 and yet we require radiation to pass through the ATR 
(and possibly OARs) to reach the PTV. 

L2 Change in the non-PTV tenns of the objective function 

A correct requirement on the OAR and ATR terms in the objective function is that they should 
be minimized to attain zero only if the weights of all the beamlets passing through these 
regions are zero. This requirement is satisfied if we use new objectivity terms of the form 



aii Dcamieis 

«>.7«= Z 2 ^fdfM- (3) 
xcAm i 

all beumlds 

Z Z wM'W. (4) 



The terms (3) and (4) will be zero only if W| is zero for all i. These terms do not allow the 
objective function to be zero through "destructive interference" effects between beamlets. 
Notice that the condition 

^_ Qltbeainlets 

X I 

is necessary and sufQcient to satisfy 

■ Zl Z ^AMj =0 with w,^0 for all i\ 
1,3 A new family of symmetry terms in ttie objective Junction 

We also introduce a family of symmetry terms to the objective function to improve the 
conformity (ICRU 1993 and ICRU 1999) of the dose distribution with positive weights for 
beamlets. The main purpose of these terms is to eliminate the interference effects that lead to 
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the appearance of negative weights. As a first example that we used for the results presented in 
this work, the term is of the form 



E • (5a) 

i 



This term complements the zero-limit for the beamlet weights in the OAR and ATR terms of 
the objective function with an equal-weight limit for all bearalets, the usual initial set of 
weights before optimization (cylindrically symmetric distribution of beams for a beam-source 
rotating around an isoccntre). The term O^^ in Eq. (5a) is positive and its minimum is zero - 
with the weights normalized to 



all bcam lea 

I 

f 



^ w^^hl — total number of beamlets , 



the minimum 0^^^ = 0 corresponds to a symmetric beamlet arrangement with w/ = 1 for all 
beams. 

A slightly different symmetry condition, less restrictive on flie beamlet weights is the 

requirement of a symmetric beam (not beamlet) distribution. Calling now wj'^ the weight of 

beamlet y inside beam i (i.e. one of the beamlets in the beam incident from the i-th gantry 
angle, and calling 

oil beamlets 
insi de be am i 

J 

the symmetry condition on the beams becomes 

oi l beam s 

I 

f 

Other possible forms of the symmetry term are 

snbeamtets 

I 

I 

all 



oi l beam s ^ . 

^^"ZK*"-!)- (5b> 



2 ^» (5c) 
# 

dlbemileis 

Z Z '^''/'W. (5d) 

xgqII conlours i 

aU btauii lcls 

x6b]| nntouis i 
■ exceptPTV 

The term (5c) is less restrictive than (5a) given fliat tibe impossibility of having all the weights 
equal to zero results in an uneven distribution of positive beamlet weights. The net effect of 
terxn (5d) is to induce an even distribution of dose deposition, everywhere regardless of contour, 
while (5e) will tend to smear the dose distribution everywhere except inside the PTV Note 
that. If desired, the contour not included in the expression (5e) can be arbitrarily chosen. The 
symmetry term included in the Qbjective function for results presented in fliis work is the even- 
beamlet term (5a). 

new objective Junction 
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0 = p^O^^+p^^O^^+p^„0^y^^;,^„0^^ (7a) 



or 



«WV t J ^toAR i 

111) beamlels ellbeamlcis 

xeAIR I I 



(7b) 



Notice that a very large value of the importance parameter p,y„ relative to the other importance 
parameters will result in an equal-weight beamlet weight distribution regardless of the shape or 
size of the contours being treated. This is a very strong condition that avoids negative weights 
but results in a non-optimized treatment plan. In practice, one needs very small relative values 
ofpjy,„ to avoid negative weights. Increasing the value of p^ym will result in a "smearing" of the 
dose distribution eliminating hot spots. Different relative values of Patr and p^y^ will result in 
different dose distributions in the regions outside the OAR and PTV. Thus, this formulation 
introduces flexibility and algorithmic control, depending on the anatomy layout and clinical 
objectives. 

2. Optimization 

2 J Condition for the minimum of the objective function 

With the modifications introduced above, the objective function remains quadratic in the 
weights.and the optimization problem for all the beam intensities is reduced to the solution of a 
system of linear algebraic equations. Given that the OAR and ATR terms of the objective 
function are identical in form, for the sake of clarity we wiU initially omit the ATR terra in the 
derivation below and reintroduce it at the end. With this in mind the objective function 
becomes 



/'aUbenmlcte V oil beamiels 



(9) 



The minimum of O is obtained by requiring 

^^di^ (10) 
Using Eq. (9) this condition results in 

/ oil bcamlets \ 
^^^^ V J :,eOAR . (11) 

Exchanging the order of summadons we can write (1 1) as: 
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' xcPTy 

ollb^ets oUbcomlett ^^^^ 



JCQOfff 

where 5« is a KrOnecker delta function; 



^ [0 otherwise * 



We define now the weight-independent arrays 

«^*».= 2: d,{x)d,(pc) (13) 

and p^<«,^^m^n, ^ ^^^^^ , ^^^^ 

where e/'^'^'^* is the dose required in region k and fi'^^''"' is a one-dimensional anay - as we 
see in expression (13), the element^ contains the sum over all points in' a region of the dose 
d/pc) deposited by beamlety! at each point in tiie region, aip'"""' is a two dimensional array - as 
wee in expression (13), the element ij contains the sum over all points in a region of the 
product of the dose di{x) deposited by bearalet i and the dose di{x) deposited by beamlet j at. 
each point x in the region in which they overlap (i.e, at each point x in the region in which both 
di{x) and d/p^ are non-*zero)» 

2,2 Solution using matrix inversion 

The condition for the minimum of O can now be written as: 

allbeamlets 

Notice that in the left side of equation (15) the weights are multiplied by the elements of a 
constant array. Reintroducing now the ATR terms we define the arrays 

«/y ^'PpTyC^r "^{PtMR^r ^PATR^r ■^Psy.)^iJ (16) 

Pj^PftvPT' ^P^'.^ (17) 
The linear system of equations to be solved becomes simply 

Y^^fit^^pj (18) 

where wsiz the (unknown) weight or intensity of beamlet >^ is a vector of coefficients that 
depends on the dose deposited by beam i within the PTV, and or^ is a matrix that describes the 
product of the doses deposited by the intersecting pairs of beamlets / and j on different organs 
(each organ with its importance coefficient). The set of optimal beam weights is obtained from 
Eq. (1) by direct inversion: 
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lg)010 



J 



(19) 



In other words, the unique solution to the large system of linear equations (18) is obtained 
quickly and accurately by inverting the array czy using standard matrix inversion routines. As 
shown later in the results section, tiie optinuzation by matrix inversion is orders of magnitude 
faster than using a numerical search algoriflun to locate a minimum. Moreover, the minimum is 
unique and global, so that local minima - a major uncertainty in numerical search techniques - 
are avoided. 

A possible complication is that system (19) may become ill-behaved if two beamlets have 
very similar dose depositions at the same spatial locations. This will result in a near-linear 
dependence of terms in the system (18) with divergent values in the matrix a"^ . The physical 
meaning is that two (or more) beamlets are nearly identical and then one or more may become 
unnecessary. This situation and the elimination of redundant beamlets are easily handled by a 
Singular Value Decomposition approach to the matrix inversion (Press et al 1992), 

23 Optimization with different importance coefficients 

The next step in planning is to choose an appropriate set of importance parameters. This is 
achieved by manual iterations of relative values of the importance parameters or through a 
programmed search, for example optimizing curves in DVH histograms. With the aid of FIDO 
either of these approaches becomes feasible given the high speed of optimization; repeat 
optimizations are not computationally onerous. In particular the "trial and error" setting of the 
importance coefficients is achieved almost interactively using this optimization routine. 

The fast re-optimization with different sets of importance coefGcients is due also to the feet 
that most of the arrays do not need to be recalculated. The arrays a*"^°"* and fi"'^""* do not 
depend on the importance coefQcients or the weights of the beamlets (they are calculated for 
unit weight beamlets), therefore they need to be calculated only once. After the first 
optimization, the arrays a*"^""* and fi"^"""' are retrieved from memory and the final equation 
arrays a and fi in equations (16) and (17) are calculated quickly by a'linear combination of 
the stored arrays and the new importance coefficients. 

3. Extensions 



3, J Region dependent importance coefficients 

In the last section we assumed that the importance coefficients for all points inside each 
contour are constants, t,e. the same for each point inside a contour. Consider now the case in 
which the importance coefficients can actually vary within the region enclosed by a contour 
(Comtz et al 2003, Shou and Xing 2004) as might be the case if tumour clonogens can be 
clearly identified with higher probability by new imaging modalities: 

Preston^ Pn^hn^QreshnS'^^ (20) 

where p is an overall constant and q{x) is a function of position. In this case, we replace fte 
definitions of equations (13) and (14) by 

Z ^r^o.M)d,ix)djix) (21) 

and =rf^<<«. 2 9r^M)d,ix). (22) 

With tiiese^ the matrix formulation of the optimization process now becomes: 
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^iS = PFrv&r APoAR^r -^PATBO^r ^P.y,n)S,j (23) 

Pj^PpTvPr -^P^v.' (24) 
The previous linear system of equations of equation (1 8) is now replaced by 

IL^Au^fij (25) 

and the optimized solution is obtained irom (25) by inversion: 

W/=IXi^i (26) 

As before, if the functional dependence of the importance coefiicients remains unaltered, a 
search of the best set of importance parameters is reduced to a search of the best set p^^^^^ in 
which case flie arrays a^^^ and fi^^^ do not need to be recalculated. Notice that the 

introduction of sub-contours is equivalent to this approach if the irnportance parameters are 
constant within a sub-contour. 

4. The circularly symmetric case 

4 J Symmetry considerations 

This simple case is remarkable because it points to the 
robustness of FIDO in a case in which the traditional 
objective function, currently in use, fails to yield 
correct results. 

Consider a target region (PTV) with a circular 
cross section inside a boundary also with a circular 
cross section as well. Both, the PTV and the exterior 
boundary, are concentric and are centered at the 
isocentre as shown in Figure 1 in which the diameters 
of the PTV and the outside boundary are 8 cm and 40 
cm respectively. Consider now an irradiation using 
arc-therapy with uniform beam intensity (not divided 
into beamlets), each beam covering exactly the PTV, 

with a constant gantry angle separation between 

adjacent beams (e.g. for 20 beams the angle separation is 18 degrees.) Given the symmetry of 
the configuration, the solution should not change if beams are rotated around the centre by an 
angle that is a multiple of the gantry angle separation between beams. As a consequence of this 
rotational symmetry, the optimal beam configuration must be such that the weights of all the 
beams must be identical . If we call W- the weight of beam i, then = Wj for all flie beams 

or, identically / wj = I for all beams, in particular for the beams with the minimum and 
maximum weights, we should find iv«^„ / w^ax = 1 - 

4.2 Failure of the conventional method 

In this case we optimized the conventional objective function 

^conventional ~ PpTV^PTV Plu^ufithmt 

with the PTV and tissue terms defined respectively in equations (2a) and (2c). 




Figure 1 
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Figures 2 and 3 plot the ratio w„i„ J w„ax for the configuration shown in Figure 1, as a 
lEunction of the ratio of the importance parameter for tissue and for PTV. A good optimization 
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should, as we have seen, yield w^ut / w^^tx = 1- Figure 2 shows the results obtained using the 
traditional objective function with a numerical search wiflxout imposing positive-weights 
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constraints. When negative weights are present then a positive*weight constrained numerical 
search yields some zero weights for which w„,i„ I = 0. If negative weights are absent, the 
results using a constrained and an unconstrained numerical search are identical. We see that the 
traditional optimization yields an incorrect result for any value oipti^zu^, the results becoming 
progressively worse as the number of beams is increased. The situation is even worse when a 
larger PTV diameter is used; as we see in Figure 4 for a PTV with a 20 cm diameter, for 60 
beams there is no value oip^ssut for which all the weights are positive! These results indicate 
that the traditional optimization yields poorer results as the overlap between beams becomes 
larger in the region outside the PTV. 

4. 3 Results using FIDO 

In this case we optimized the FIDO objective function 



with the PTV and tissue terms defined respectively in equations (2a) and (4). No symmetry 
term, of the type of equations (5), was included for this calculation. As shown in Figures 3 and 
5, when the FIDO method is used, once acquires a minimally significant value {passue / 
Pptt^ 0-01) for tiie tissue term to become relevant in the objective function (i.e. when FIDO 
"kicks in"), all the beams have the same optimized weight within the precision of the 
calculation. No difference is found in the case of a PTV with a larger radius. 
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These results are very important because they show that FIDO is not only much faster and 
accurate than traditional optimizations but also that there are limiting test cases in which FIDO 
yields correct results whUe the conventional objective Junction (Le. the standard optimization 
currently in use) provides erroneous solutions. 

5« Optimization Results 

Below we present several sets of preliminary results in two dimensions (2D) obtained for a 
head and neck case, a prostate case and an "interlocking rectangles" phantom. In each case, we 
include calculation and optimization times for primary KERMA calculations as they were 
obtained on a single-processor desktop PC. The program was written in C# using the Microsoft 
.NET environment. No effort has been devoted to maximize the speed of the matrix inversion 
procedure. 

In all cases the number of gantry angles used is evenly distributed over a full 360 degree 
circle around the isocentre. In turn, each beam is evenly divided into beamlets of the specified 
width resulting in the total number of beamlets quoted. The source to axis distance (SAD) is 
100 cm. In each case we also present the cumulative dose-volume histogram (DVH), as well as 
a set of iso-KERMA contours. 



5.1 Head and neck case 



In Figures 6 and 7 we present the isodose curves and the dose-volume histogram (DVH) 
for a neck tumour partially wrapped around the spinal cord with a linear dimension of 
approximately 11 cm. The target was iiradiated with 60 equally spaced beams divided into 
beamlets 2 mm wide, yielding a total of 2968 beamlets. The FIDO optimization of the 2968 
independent weights using a Matlab matrix inversion routine takes 5.13 seconds on a single 
processor desktop PC. The resulting dose distribution within the PTV is homogeneous within 
4>A 2% while 95% of the OAR receives a dose that is less than 5% the dose in the PTV. 



Figure 6 
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5.2 Interlocked square phantom case 

In Figures 8 and 9 we present the isodose curves and the dose-volume histogram (DVH) 
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Figure 9 
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for a challenging "interlock" phantom. The PTV, with a linear dimension of approximately 
25 cm is wrapped aroimd half of a rectangular OAR. The target was irradiated with 40 equally 
spaced beams divided into beamlets 3.75 mm wide, yielding a total of 2326 beamlets» The 
FIDO optimization of flie 2326 independent weights using a Matlab matrix inversion routine 
takes 2.53 seconds on a single processor desktop PC. The resulting dose distribution within the 
PTV is homogeneous within +A3% while 95% of the OAR receives a dose that is less than 
20% the dose in the PTV. Notice that the 97% isodose is very highly conformal and that any 
iso-Kenna for values greater or equal to 103% is not observed; notice as well the sharp 
gradient in the region between the PTV and the OAR. 

In both cases above, i.e. the optimization of 2,968 beamlets and 2,326 beamlets, a Matlab 
minimization search routine with gradient input, using the standard objective lunction and 
constrained to positive weights, nmning under identical conditions as in the cases above, could 
not find a global minimum in less than 4 hours (at that time die program was shut down). 

5.3 Prostate with overlapping rectum case 

In Figures 10, 11 and 12 we present the isodose curves and the dose-volume histogram 

(pVH) for a prostate PTV that overlaps 
the rectum. The PTV, with a linear 
dimension of approximately 7 cm is 
wrapped around half of a rectangular 
OAR. The target was hradiated with 60 
equally spaced beams divided into 
beamlets 4 mm wide, yielding a total of 
1002 beamlets. The FIDO optimization 
of the 1002 independent weights using 

a Matlab matrix inversion routine takes 

0.24 seconds on a single processor 
desktop PC. In Figure 10 we present a panoramic view of the slice including the femoral heads. 
In Figure 11, we present a close up view of the PTV-rectum region, and in Figure 12 we 
present the Dose- Volume Histogram - notice that in the case of the OAR we plot the DVH for 
the portion of the OAR outside the DVH (i.e. the portion that needs to be shielded). The 




Figure 10 
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Figure 12 



resulting dose distribution within the PTV is homogeneous within +4% and -3% with the white 
dots in side the PTV m Figure 11 corresponding to 103% isodoses. Notice again the sharp 
conformality to the PTV of the 97% iso-KERMA as well as the sharp gradient as one enters 
the OAR portion outside the PTV. 

. In this case, i.e. the optimization of 1,002 beamlets, a Matlab minimization search routine 
with gradient input, using the standard objective function and constrained to positive weights, 
running under identical conditions as in the cases above, took 3 hours and 27 minutes to 
converge to the global minimum. 

The results above display the power of FIDO to perform optimizations extremely quickly 
and to provide very conformal dose distributions - both important towards the goal of adaptive 
radiotherapy. 
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6. Conclusions and future work 

We have developed a fast and robust technique to find a global minimum that yields excellent 
results for the inverse Optimization problem for the radiation treatment of tumours with IMRT, 
using large sets of non-negative intensity-modulated beamlets. 

Work is currently in progress on a full implementation FIDO in our treatment planning 
system that uses collapsed cone convolution method for dose calculation (Ahnesjo 1989). 
Work is proceeding as well on an interactive 3D iniplementation that will permit "dose- 
adaptive radiotherapy" using CT images (Mackie et al 1999, JafFray et al, 2002): The dose 
distribution can thus be reconstructed for each treatment fraction to converge to the planned 
optimization goals. 
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FIGURE CAPTIONS 

Figure 1. Concentric circular contours for the circularly synunetric case. 

Figure 2. Standard optimization results for the weight ratio w^Jw^ as a function of 
importance parameters for a circular PTV with a 4 cm radius. 

Figure 3. FIDO results for the weight ratio w^Jw^^ as a function of importance 
parameters for a circular PTV with a 4 cm radius. 

Figure 4. Standard optimization results for the weight ratio w^ijw^^ as a function of 
in^ortance parameters for a circular PTV with a 1 0 cm radius. 

Figure 5, FIDO results for the weight ratio w^iJw^ as a function of importance 
parameters for a circular PTV with a 10 cm radius. 

Figure 6. Isodose contours for a head and neck case. 

Figure 7. Dose-volume histograms PVH) for the head and neck case of Figure 6. 

Figure 8. Isodose contours for an interlocked phantom. 

Figure 9. Dose-volume histograms (DVH) for the interlocked phantom of Figure 6. 

Figure 10. Isodose contours for a prostate case - panoramic view. 

Figure 11. Isodose contours for a prostate case - close-up view. 

Figure 12. Dose-volume histograms (DVH) for the prostate case of Figures 10 and 1 1. 
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